Transition to Chaos in the Self-Excited 
System with a Cubic Double Well Potential 
and Parametric Forcing 



Grzegorz Litak, 81,1 Marek Borowiec, a Arkadiusz Syta, b 
Kazimierz Szabelski a 

1 Department of Applied Mechanics, Technical University of Lublin, Nadbystrzycka 

36, PL-20-618 Lublin, Poland 

b Department of Applied Mathematics, Technical University of Lublin, 
Nadbystrzycka 36, PL-20-618 Lublin, Poland 



Abstract 

We examine the Melnikov criterion for a global homoclinic bifurcation and a possi- 
ble transition to chaos in case of a single degree of freedom nonlinear oscillator with 
a symmetric double well nonlinear potential. The system was subjected simultane- 
ously to parametric periodic forcing and self excitation via negative damping term. 
Detailed numerical studies confirm the analytical predictions and show that tran- 
sitions from regular to chaotic types of motion are often associated with increasing 
the energy of an oscillator and its escape from a single well. 



1 Introduction 



A system of nonlinear stiffness having a square displacement force term, and 
simultaneously, excited externally or parametrically has been a subject of stud- 
ies for many years [1,2,3,4,5,6,7,8,9,10]. It has found numerous applications 
in mechanical engineering [1,2,3,4,5,6,7,8] and control theory [9,10,11,12]. It 
was also one of intriguing examples of simple nonlinear systems which showed 
complex behaviour including chaotic oscillations [13]. Another important, par- 
tially separated, class of systems can be determined on the basis of nonlinear 
damping. Such damping can, in some systems, change the sign depending on 
velocity or displacement values, and provide excitation energy to the examined 
system. These, so called, self-exited damping terms are often used to describe 
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systems with dry friction, bearings lubricated by thin layer of oil, shimming 
in vehicle wheels or chatter in cutting process, [1,2,3,14,15]. More recently 
such a nonlinear damping force has also been considered [16] in modelling of a 
modern vehicle suspension system due to electro- or magneto-rheological fluid 
damping where it is causing a hysteretic effect. In this model [16] the authors 
used a self-excited term of Rayleigh type and Duffing double well potential. 

This class of a system has also been investigated by Siewe at al. [17,18] in 
their works on $ 6 -Van der Pol oscillators. In the first paper [17] they applied 
Melnikov approach and performed numerical simulations in case of a poten- 
tial based on polynomials of even orders (4th and 6th orders: $ 4 and $ 6 type 
potentials respectively) and an external excitation. On the other hand in sec- 
ond paper [18] the triple well potential based on polynomial of the 6th order 
was studied in details for a dynamical system with external and parametric 
excitations. The Melnikov approach has been also considered in systems with 
a single well potential and various nonlinear damping terms. Among a large 
number of papers conducted that research we report [19,20,21]. Litak et al. 
[19] included the Rayleigh self-excited term to describe dynamics of the Froude 
pendulum, Trueba et al. [20] considered nonlinear damping following a power 
law in velocity (a(x) n , n = 2,3) for various systems, and finally Awrejcewicz 
and Holicke [21] examined the effect of dry friction in a stick-slip system with 
the Duffing potential. 
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Fig. 1. The double well potential V(x) = 5x 2 /2 + ^y\x\x 2 /3 for three values of 7 
(7 = 0.8, 0.5, and 0.2) and 6 = -1.0 . 

Note, that all above systems can be still regarded as simple ones but its com- 
bined nonlinear answer on excitation is complicated and deserves detailed 
investigations. 

In this note we shall examine transition from regular oscillations to chaos 
in a simple, one degree of freedom, system subjected to parametric and self 
excitations with a square albeit symmetric stiffness: 

x + a(—l + x 2 )x + (5 — ficos2cot)x + j\x\x = 0, (1) 
where x is a displacement a(— 1 + x 2 )x is nonlinear damping, —fixcos2uit is 
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a parametric excitation while ^\x\x and Sx are square and linear force terms. 
It is worth mentioning that, the corresponding potential can be described by 
two polynomials of the 3rd order (<3> 3 type) for x > and x < 0, respectively: 



5x 2 

V(x) = + 
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(2) 



and plotted in Fig. 1 for 7 = 0.2, 7 = 0.5 and 7 = 0.8. In spite of appearing 
|x| in Eq. 2 the function V(x) is of C 2 class at x = because it is approaching 
to 0, for x — > 0, as ±x 3 . 

The present paper consists of three sections and two appendixes. After a short 
introduction in the present section (Sec. 1) we will examine the above system 
analytically and numerically in Sec. 2. There we will especially focus on a tran- 
sition and a corresponding scenario from regular to chaotic vibrations. In Sec. 
3 we will discuss the numerical simulations verifying the analytical findings. In 
the last section (Sec. 4) we will provide the summary and conclusions. In the 
appendixes we will present different treatment of a self-excitation term and 
its consequences to the applied analytical methods (Appendix A) and show 
the Melnikov integration procedures in details (Appendix B). 



2 Melnikov Analysis 

We are starting our study from the second order equation of motion (Eq. 
1). After transforming it into two differential equations of the first order, 
a standard procedure in examining homoclinic transition will be applied to 
look for stable and unstable manifolds and their possible cross-sections in 
presence of weak excitations and damping. Therefore we have introduced a 
small parameter e to the above equations enabling these terms to be switched 



where ae = a and jle = [i. Note, that this is not a unique way of splitting the 
initial differential equation of the second order (Eq. 1) into two equations of 
the first order. 

The other way is connected with different treatment of the self-excitation 
term. Basing on 'fast' and 'slow' variables (x,w) identification in Van der Pol 
fashion [24] we can write: 



on [4,22,23]: 



x — v 



(3) 



v = — 5x — j\x\x + e a(l — x 2 )v + xp, cos (2ut) , 
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Fig. 2. Comparison of v and w time histories for 7 = 0.2 (Fig. 2a) and 7 = 0.8 (Fig. 
2b). Other system parameters: a = 0.1 uj = 0.45,5 = — 1.0, // = 0.6. (Fig. 2c) Phase 
portraits (in x-v space) of for both cases (7 = 0.2 and 0.8). 



w = x + ea Ix — — 



(4) 



This possibility has been examined thoroughly in Appendix A. Note, the dif- 
ference between w(t) and v(t) depends on the influence of the self-excitation 
term leading to relaxation oscillations in the Van der Pol system. The time 
histories of w and v (Figs. 2a-b) show that this term could more influential 
for smaller 7 (7 = 0.2, Fig. 2a) while a case with larger 7 (7 = 0.8, Fig. lb) is 
apparently negligible. For better clarity we plotted phase portraits (Fig. 2c) 
for both cases showing that they correspond to steady state vibrations located 
in different regions of the phase plane. This is associated with differences in a 
potential shape (Fig. 1). 

In the first case vibrations have been realized around the minimum of the 
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potential V(x) (x ~ 5.0) inside one of its wells, while in the second one 
vibrations are located around x = between x = —2 and 2. As in the examined 
example the self-excitation Van der Pol term: 



VdP(x,x) = a(-l +x 2 )x 



(5) 



is changing its sign at x — ±1.0. One can easily see that for 7 = 0.2 the system 
behaves if it possessed renormalised nonlinear damping term which does not 
change its sign. Moreover this term is positively defined. Thus the system, due 
to the shape of the external potential (Fig. 1), does not have any relaxation 
character typical for simple Van der Pol oscillator and the distinction of 'slow' 
and 'fast' variables is not relevant here. This fact has its reminiscence in a 
large difference in v(t) and w(t) time histories in Fig. 2a. The transformation 
v(t) — > w(t) can be better argued in the the second examined case (7 = 0.8). 
There the relaxation vibrations are present but in this case the effect of self- 
excitation turns out be small (Fig. 2b). 

Note, that the unperturbed Hamiltonian H° is the same for both cases (with 
v or w variables) and reads: 



The potential function V(x) (Fig. 1) has the local peak at x = 0. The exis- 
tence of this point with a horizontal tangent makes global homoclinic bifurca- 
tions, including transition from regular to chaotic solution, which may occur 
in the system. Obviously, at this point the system velocity reaches zero value 
(v = 0) (Fig. 1). Thus, according to our potential gauge, the total energy has 
only its potential part which is also zero: 



Transforming Eqs. 4,6 for a constant energy, chosen here as zero, (Eq. 7) we 
obtain the following expression for velocity: 



H° = - + V(x). 



(6) 



E = V(x = 0) = 0. 



(7) 




(8) 



from which 




(9) 
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where to represents an integration constant. As a result of integration (Eq. 9) 
we get so called homoclinic orbits (Fig. 3) parametrised by time t: 



35/ ^ 2 (v=s(t-toy 



x * = x *(t- t ) = ±ri M - tanh 2 . 



v* =v*(t-t ) = T 



3S ^S tanh(^p>) 
2 7 cosh 2 (v^z!o))' 



(10) 



where '+' and '— ' signs are related to I and II orbits (Fig. 3), respectively. 




Fig. 3. Homoclinic orbits (I and II) of unperturbed potential V(x) = —x 2 /2+\x\x 2 /3 
(for 5 = — 1 and 7 = 1) in a phase plane (x,x). For t — > ±00 we get (x,x) — > 
(x ,xo) = (0,0). 



Note, the central saddle point xq = is reached in time t corresponding to 
+00 and —00, respectively. 

In case of perturbed orbits W s (a stable manifold) and W u (an unstable 
manifold) the distance between them is given by the Melnikov function M(t ) 
[22,23,26]: 



+00 

M(t„) = J h(x*,v*) A g(x*,v*)dt 



(11) 



where the corresponding differential forms h is the gradient of unperturbed 
Hamiltonian (Eq. 6): 



h = (8x + 7|a;*|x*) da; + vdv, 



;i2) 



and g is a perturbation form (Eq. 3) to the same Hamiltonian (Eq. 6): 

g = {jxx cos 2uot + a (l — x 2 ^j v^j dx (13) 
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Both differential forms are defined on homoclinic manifold (x,v) = (x*, 
(Eq. 10, Fig. 2). 
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Fig. 4. Critical value of [i c , for a = 0.1, versus u for 6 = —1 and various 7 (7 = 
(Fig. 4a), 7 = 0.5 (Fig. 4b), 7 = 0.2 (Fig. 4c)) and versus 7 (Fig. 4d). 



Naturally, from the above (Eqs. 11-13) the Melnikov integral is given by: 
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Fig. 5. Phase diagrams and Poincare sections (simulations were based on Eq. 1) 
chosen sets system parameters a = 0.1, 5 = —1.0 and u, 7 (corresponding 
values are indicated in figures) denoted by points in Figs. 4a-d. The corresponding 
maximal Lyapunov exponent is (Fig. 5a): Ai = —0.0279 for '1' and Ai = —0.0238 
for '2', (Fig. 5b): Ai = 0.0903, (Fig. 5c): Ai = -0.0978), (Fig. 5d): Ai = 0.0812, 
(Fig. 5e): Ai = -0.1513, (Fig. 5f): Ai = 0.1202, (Fig. 5g): Ai = -0.1974, (Fig. 5h): 
Ai = 0.0529, respectively. 



+00 



M(t ) = / v* (t- t ) {jlx* (t-t ) cos (2ut) + a 1 — x* (t-t Q )]v*(t-t )} dt 



—00 
+00 



v*(t) {jlx*(t) cos (2u(t + t )) + a [l - x* 2 {t) 



v*(t)\ dt. 



(14) 



After substituting x*(t) and v*(t) given in Eq. 10 we obtain: 
M(*„) = fih{t ) + al 2 , 



(15) 
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Fig. 5. Continuation. 

where the component integrals I\(t ) and J 2 are defined 
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I^to) = / w*(t)x*(t)cos(2cj(t + t ))dt 



(16) 



and 



(l-x* 2 (t))v* 2 (t)dt. 



(17) 



After evaluation of the above integrals (see Appendix B: Eqs. B.17 and B.4 ) 
we get the exact forms of h(to), and 



T , . 12tt5 2 J 1 (1 + 4cu 2 ) . 

Ji(*o) = — 2 r T7o' — r sm ( 2 ^o) 

7^ smh(2o;7rJ 



(18) 
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Fig. 6. Stable and unstable manifolds in presence of a perturbation component 
(e = 0.05). Parameters used in calculations (Eq. 2): uj = 0.4, 7 = 0.2, 5 = —1.0, 
a = 0.1 and various excitation amplitude j2 = 0.87 (Fig. 6a), 0.94 (Fig. 6b), 1.8 
(Fig. 6c), 2.7 (Fig. 6d). Note that \i = eft, a = ea (see Eq. 3). Arrows in Fig. 6c 
and d indicate the crossing points. 
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Surprisingly the other chose of system variables (x, w) lead to the same form 
of Melnikov function (Appendix A), Eqs. 14-17. That means that both per- 
turbation procedures give the same result in the first approximation. This is 
partially due to the fact that in both cases we started from the same set of 
equations of unperturbed Hamiltonian H° (Eqs. 2,6). 



Finally, the condition for a transition to chaotic motion as a global homoclinic 
transition corresponding to a horse-shoe type of stable and unstable manifolds 
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Fig. 7. Time histories x(t) obtained by simulation of Eq. 2 u = 0.4, 7 = 0.2, 
5 = -1.0, a = 0.1, [i = 0.94 and three values of e = 0.05 (Fig. 7a), e = 0.50 (Fig. 
7b), e = 1.00 (Fig. 7c). Note that \i = e/2, a = sa (see Eq. 3). 



cross-section, can be written as 



V M(f ) = and ^ 0. (20) 

t Mo 



The above integrals (Eq. 18,19) together with the last condition (Eq. 20) yields 
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Fig. 8. Maximal Lyapunov exponent Ai (Fig. 8a) and the minimum point of dis- 
placement point Xmi n (Fig. 8b) of steady state vibrations. Figs. 8d-e correspond to 
bifurcation diagrams for 7 = 0.8, 0.5 and 0.2 respectively for a = 0.1, 5 = —1 and 
<jj = 0.4. Simulations were based on Eq. 1. 

a critical value of excitation amplitude \x c : 
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for which stable and unstable manifolds cross. 

Namely we get the critical amplitude \i c versus frequency u>, and parameter 
7 which is plotted in Figs. 4a-d. Above this value /i > /i c the system transit 
through a global homoclinic bifurcation which is a necessary condition for ap- 
pearance of chaotic vibrations. Interestingly, larger 7 leads to smaller // c (Fig. 
4d) Note also, that the condition (Eq. 21) is based on results of perturbation 
procedure in its lowest order approximation. Equation 21 is the main results 
of our investigation. 



3 Numerical simulations 



To verify the analytical findings we have performed series of numerical sim- 
ulations of Eq. 1. In Fig. 5 we show phase diagrams and simultaneously 
corresponding Poincare sections stroboscopic points collected with frequency 
il = 2uo for chosen sets initial system (Eq. 1) parameters: uo = 0.4, 7 = 0.8 
H = 0.02, 0.10 (Fig. 5a) and n = 0.116 (Fig. 5b); cu = 0.4, 7 = 0.5 and 
H = 0.160 (Fig. 5c) and n = 0.185 cu = 0.4, 7 = 0.2 (Fig. 5d); n = 0.87 
(Fig. 5e); and /j, = 0.94 Fig. 5f); u = 0.45, 7 = 0.2 and n = 0.90 (Fig. 5g); 
and fi = 0.97 Fig. 5h). The structure of the examined attractors as well as 
calculated Lyapunov exponents enables to classify the dynamics of the sys- 
tem (see Fig. 5 and figure caption). For comparison with analytical results 
we indicated the simulated cases by points, showing the types of vibrations: 
R (regular) and C (chaotic). For the last case for relatively small 7 (7 = 0.2 
in Fig. 4c, and Figs. 5 e-h) one can see some discrepancy between simulated 
data and analytical results. Analytical results indicate that the homoclinic 
transition take place for /i c > 1 but simultaneously for the same system pa- 
rameters numerical results states that fi c < 1 is enough to transit into chaotic 
vibrations (/i c ~ 1.5, Fig. 4c). This could be effect of relatively large fi in 
perturbation procedure (Eq. 3). In the above calculations we used the same 
initial conditions (x in ,v in ) = (0.45,0.1). 

The interesting discrepancy between analytical (Fig. 4c) and numerical simu- 
lation (Fig. 5e-f) results deserved some additional comments. In this aim we 
plotted, in Fig. 6) our results on simulations of perturbed system (Eq. 3) re- 
lated to stable and unstable orbits. Because the Melnikov theory formulated 
in the lowest order approximation is valid in the limit of small e (e — > 0) we 
decided to use a relatively small value of this parameter e = 0.05. Other pa- 
rameters used in calculations (Eq. 3) are as follows uo = 0.4, 7 = 0.2, 8 = —1.0, 
a = 0.1 and various excitation amplitude p, = 0.87 (Fig. 6a), 0.94 (Fig. 6b), 
1.8 (Fig. 6c), 2.7 (Fig. 6d). Note, in cases of Fig. 6a and b, which correspond 
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according to Figs. 4c, 5e, 5f, to regular and chaotic motions for e — 1 (original 
equation Eq. 1). Obviously there are no crossing points between the stable 
W s and unstable W u manifolds. After comparing these results it is obvious 
that in this particular case the final result of our consideration could depend 
on e. In fact the difference appears between the limit e — > (or e = 0.05 in 
Figs. 6a,b) and e = 1.0, what is not surprise taking into account relatively 
large value of /t (jl ~ 1.0 Fig. 4c). Interestingly for the present value of e ( 
e = 0.05) our simulations (Fig. 6c and d) confirm analytical predictions (see 
the solid curve in Fig. 4c). As /i > /j, c 1.5 in both cases we observe crossing 
points of the stable and unstable orbits (indicated by arrows in Fig. c and 
d). Evidently increasing jl causes faster meeting of corresponding curves W s 
and W u , as one can see Fig. 6c and Fig. 6d differ by winding numbers N of 
motion necessary for such crossing: N — 1 for Fig. 6c and N = 2 for Fig. 
6d. We have also done additional tests for other system parameters including 
different 7 = 0.5, 0.8 (as in Fig. 4a-b and Figs. 5a-d) and uj = 0.45 (as in 
Fig. 4c and Figs. 5g-h). In all cases we have also got agreement with the cor- 
responding analytical curves (Fig. 4). In the next figures (Fig. 7a-c) we show 
the evolution of the oscillation with increasing e. Here we have plotted time 
histories for different values of e. Starting from vicinity of the saddle point 
we obtained regular history (Fig. 7a) for e = 0.05. For larger value of control 
parameter e (e — 0.5 in Fig. 7b) we observe enlargement of the amplitude of 
these motion and finally its transition to chaotic motion for a large enough e 
(e = 1.0 in Fig. 7c). 

For better clarity we decided to show the transition to chaos through the 
maximal Lyapunov exponent Ai and bifurcation diagrams for three values of 
7 (7 = 0.2, 0.5, 0.8) versus \x (Fig. 8). Positive value of Ai in Fig. 8a( see the 
region ~ 0.15 < /i <~ 0.45 for 7 = 0.8) as well as black regions in correspond- 
ing bifurcation diagrams (Fig. 8c-e) imply transition to chaos. This is a way 
of numerical identification of critical value /i c regarding a local bifurcation. In 
case of 7 = 0.8 and 0.5 we observe the confirmation of previous findings for nu- 
merical (Poincare maps Figs. 5a-d) and analytical results (Fig. 4a-b) while the 
case of 7 = 0.2 is different. Here again we observe some disagreement between 
numerical (Figs. 5e-h) and analytical results (Fig. 4c). Lyapunov exponent 
method gives smaller value of /i c . The results, obtained while simulating Eq. 
1, show the actual critical values of [i for transition to chaotic vibrations. Note 
also, that in many cases the observe the chaotic oscillations existing in both 
potential wells (with positive and negative x). However this behaviour is not 
a rule because Figs. 5g,d show the transition to chaotic behaviour inside the 
single well x > 0. To explore this effect more systematically we plot (in Fig. 
8b) x min versus \i and three values of 7 = 0.2, 0.5, 0.8. It appeared that for 
7 = 0.8 escapes and returns of the considered system to a single well oscillatory 
motion can happen for an interval of \x parameter which is close but slightly 
below the value transition to chaotic vibration. This means that the system 
has already possessed double well attractor just before its transition to chaotic 
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behaviour. Note, in our numerical calculations, to have some direct insight into 
the system capability of transition to chaotic vibration and its escape from the 
a single well and to minimize the time of our calculations, the initial condi- 
tions at the beginning (// = 0) were assumed to be (x in ,v in ) = (0.45,0.1) and 
continuing the calculations we used final displacement and velocity obtained 
for the former value of fi (xfi,Vfi) as initial conditions for any new value of 
[i. We also note, that this system resembles the other examined by Thompson 
[13] with the same type of nonlinearity of square type but without Van der 
Pol term and external forcing instead of our parametric one. In his study he 
has got the condition for a global homoclinic transition just before the escape 
from the single potential well. In our case the effect of self excitation, con- 
nected with the Van der Pol damping term is stronger for 7 = 0.8 (Fig. 1). 
This effect can be deduced from the excitation term itself for x — > we have 
VdP(x,x) — > — ax which means negative damping if the system moves close 
to x = 0. The extra energy generated by such a negative damping makes it 
possible to overcome the potential barrier between two symmetric potential 
wells much easier comparing the case without self excitation effect. On the 
other hand, the size of attractor for given 7 (Fig. 5) is solely determined by 
the shape of the nonlinear potential (Eq. 2, Fig. 1). For 7 = 0.8 we have also 
found that for /i > /i c a transition chaotic vibrations was preceded by transient 
chaotic behaviour. The steady state motion undergo escape from a potential 
well with formation of symmetric attractor spreading on two potential wells. 
With increasing \x the symmetry of attractor was broken and the system tran- 
sit a series of period doubling bifurcations as was discussed in earlier papers 
[7,8,26]. 



4 Summary and Conclusions 

In summary we have studied conditions of a global homoclinic bifurcation in 
a double well potential -van der Pol system with parametric excitation. Such 
a bifurcation correspond to transit chaotic behaviour of the system and with 
some further increasing of the excitation amplitude can lead to the permanent 
chaos [7,8,14,26]. 

Using the Melnikov method we have got the analytical formula for transition 
to chaos in a one degree of freedom, system subjected to parametric excitation 
with a non-symmetric stiffness with self-excitation term. In our case this effect 
is mutually introduced through the Van der Pol damping and parametric ex- 
citation terms. Our analytical results are consistent with direct computations 
on homoclinic orbits. 

Note that our vector field (Eq. 3) is of C 1 class due to non-continuity of the 
second derivative at x = line or a piece-wise C 2 smooth system. In fact 
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he standard theory [22] assumes that expansion in Taylor series in e is good 
to second order at least (for uniform bounds on e 2 term). This requires that 
h(x, x) and g(x, x) are of C 2 class with respect to x and x, and the Hamiltonian 
is of C 3 class. However for specific cases the theory can be applied for weaker 
assumptions. For any non-smooth systems one should check that piece-wise 
smooth solutions can be assembled and indicate any jumps in derivatives that 
can occur [25]. 

Let us write the forms h(x,x) and g(x, x) (Eqs. 12 and 13) as vectors: 



Now, the Mielnikov function M(to) can be treated as a projection of the vector 
= [— hi-, hi] into the g direction (a scalar product), Namely 



Note that the vector element h£ (h^ = hi) is projected out the Melnikov 
integral. Moreover the same argument applies to any x derivative of hi making 
its non-smooth behaviour at x = unimportant for the Melnikov theory 
application. In other words this non-smoothness at x = is not likely to 
produce any jump to the homoclinic orbit (Fig. 3). 

Interestingly the families of function fj, c (u) plotted against uj scales as ~ 7~ 2 
leading to small /i c for relatively large 7 (7 = 0.8) and much larger /i c for 
7 = 0.2. To confirm these results we have performed numerical simulations 
showing corresponding phase diagrams, Poincare maps Lyapunov exponents 
and bifurcation maps. The Lyapunov exponent has been calculated using the 
algorithm provided by Wolf et al. [27]. 

We have noticed some discrepancy (due to relatively large values of [i c : \\i c j S\ > 
1), between simulated data and analytical results, especially, in case of small 
7. It seams that our Melnikov analysis provides results which are correct up 
to the first order. The results can be possibly improved in higher order ap- 
proximations [28] The other possibility is the stronger influence of the Van 
der Pol term on the system dynamics. However, keeping the Van der Pol term 
unchanged, we have not analyzed this effect so deep as it deserves leaving this 
aspect to future studies. 



h= [(5x + -f\x*\x*) ,v] = [hi,h 2 ], 




(22) 




(23) 



—00 



16 



Acknowledgements 

This paper has been partially supported by the Polish Ministry of Science 
and Informatization. GL would like to thank the Max Planck Institute for the 
Physics of Complex Systems for hospitality. Authors would like to thank Prof. 
P. Holmes for discussion and unknown reviewers for valuable comments. 

Appendix A 

The starting equation of motion (Eq. 1) is given as a differential equation of 
the first order. To perform further analysis of the system it must be split into 
two equations of the second order. The usual way based on time derivatives 
of displacement and velocity (x,v) was discussed in Sec. 2 (Eq. 2). The other 
concept (Eq. 3), connected to division on 'slow' and 'fast' variables leads to a 
another pair of equations: 



w = —5x — j\x\x + ejl (cos (2ut)) . 

The above splitting does not effect the unperturbed Hamiltonian H° which 
is of the same form as for (Eq. 2) Thus w has the same meaning of v but 
perturbations are now appearing in both equations instead of one. 

Now, the gradient of unperturbed Hamiltonian: 

h = (Sx + 7|x* \x*) dx + vdv, (A. 2) 

while g' is a perturbation form 



defined on corresponding stable or unstable manifolds (x,v) = (x*,v*) (Eq. 
10, Fig. 2). The Melnikov function: 




(Al) 




(A.3) 



+ 0O 




(A.4) 



— oo 



+ 0O 



*3 



) 



) 



fix*v* cos (2u(t + t )) + a (6x* + 7a;* 2 ) I x 



X 



dt 
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can be evaluated after substituting x*(t) and v*(t) by formulae given in Eq. 
10: 

M'(t )=fiI' 1 (to) + aI^ (A.5) 
where 

oo 

I[(t )= j dt v*(t)x*(t)cos(2u(t + t )) (A.6) 

— oo 

and 

I' 2 =] dt (6x*(t) + 1X *\t)) (x*(t) - ^) . (A.7) 

— oo * ' 

Note the expression for /((to) (Eq. A.6) coincides exactly with Ii(t ) analyzed 
in Sec. 2 (Eq. 16) while I' 2 can be transformed exactly to I 2 if one make 
use of the definition of the homoclinic orbit parametrisation (Eqs. 7-8) and 
substituting: 

dv 

— = —5x — 7|x|x (A. 8) 

into the expression Eq. A.7 and integrating it by parts: 

—00 v ' 

00 

= J dt(l- x* 2 (t)) v* 2 (t) = I 2 . (A.10) 
—00 

Thus the above Melnikov function form M'(to) coincides exactly the the form 
of M(t ) obtained Sec. 2 (Eqs. 14-19). 

Appendix B 

Evaluation of the integral I2 is straightforward. After substitution x*(t), and 
v*(t) (Eq. 10) we have: 

00 

I 2 = J (l-x* 2 (t))v* 2 (t)dt (B.l) 
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Fig. B.l. Deformed contour integration schema and imaginary poles. 



9^3 J tanh 2 



-5f 



47 •/ cosh . , 

— 00 \ 2 



9 5_ 2 
4 7 2 \ v ' 



1 - -— I tanh 2 



-St 



1 dt. 



and simple algebraic manipulations: 
2r 



t 



tanh r = £ 



we obtain 



"2 7 2 



fB.2) 



/e 2 (i^ 2 )(i-^(^-i) 2 )^. (b.3) 



Finally the result of integration one has a following expression: 
8 2 V^S ( 6 72<5 2 \ 



7' 



5 707^ 



On the the hand the integral I\ can be written as follows 



(B.4) 



J 1 (t )= / w*(t)x*(t)cos(2w(t + to))di 

—00 

9 5 2 J t tanh v 2 

= ~2t^ 



2 cos h 2 (^) 



1 — tanh 



-5^ 



(B.5) 

cos(2cu(t + t ))dt. 



what can be expressed further in terms of more elementary integrals Jj 4 and 
If 



h(t ) 



35 2 



2 7 2 

which are given by: 



(If - 3if ) sin(2cjt c 



(B.6) 
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. °r tanhr / Acur \ , 

— oo N ' 

tR °f tanh 3 r . / Aur \ , 



(B.7) 
(B.8) 



Note that the second integral If can be easily obtained from the first one Jj 4 
through integration by parts 

oo 12 ( 8 



(B.9) 



where 



a; 



4a; 



(B.10) 



While Ii should be calculated using the residue theorem 

r N 
ff(z)dz = 2iriJ2Res[f(z),z k ], 

J i 1 



k=i 



(B.ll) 



where 



lm— 1 



Res[f(z),z k ] = j—^ hm _ [(z - z k Tf{z)\ 



(B.12) 



In our case 



„. , exp(z) - exp(-z) iiuz 

./(•?) = 4 , : : ex P ■ 



(exp(z) + exp(-z)) 3 y/^S 



(B.13) 



where on the real axis (Fig. B.l) Kez = r: 

T „. . tanhr . / Aujt \ 
lm/ [z) = r^5— sin 



cosh r 



-5 



(B.14) 



The multiplicity of each pole of the complex function f(z) (Eq. B.13): 



:,, = [^ + nk)\ foi /, = 1.2.3. 



(B.15) 
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can be easily determined as m = 3. After summation of all poles (Fig. B.l) 
we get: 



o sinh 



The result of the above analysis can be written in a compact way: 

/l(to) ^^l±^) 8in(2 ^). (B ,7) 
7 Z smh(2u;7rj 
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